import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.optimize import linprog

# 读取文件
df = pd.read_csv("E:/data.csv")

# 本程序运用Python内置的linprog函数对问题进行求解，函数声明为：
# linprog(c, A_ub, b_ub, A_eq, b_eq, bounds, method)
# c——目标函数的决策变量对应的系数向量（不限定是行或列向量，下同）
# A_ub, b_ub——不等式约束条件的变量系数矩阵和常数项矩阵（必须是<=形式）
# A_eq, b_eq——等式约束条件的系数矩阵和常数项矩阵
# bounds——表示决策变量定义域的n*2矩阵，None表示无穷。默认为（0，None）。另：Python中常用float('-inf')表示负无穷
# method——调用的求解方法，默认为highs
# 返回值——有多个成员，常用成员有：x为最优情况下各决策变量应有的取值（最优解），fun为函数最小值，nit为迭代次数。调用方式举例：
# result = linprog()
# x = result.x（似乎就像结构体一样）
risk_thresholds = np.array([])
earnings = np.array([])  # 存放a与相应最值
earnings_pure = np.array([])  # 存放净收益最值
solutions = []  # 存放各结果对应最优解
c = np.array([0.05])  # 额外再加上银行利率，没有相应的费率需要减
c = np.append(c, (df['r_i']-df['p_i']) * 0.01)
# 1.linprog函数只能求解最小值问题。如果希望求最大值（如本题），可以将所有系数变为相反数
c = -1.0 * c
M = 10000  # 假设出的M的值
A_ub = np.diag(np.append(np.array([0.0]), df['q_i'] * 0.01))  # 每个不等式约束中，有且仅有一个变量系数不为0，
# 另外，对于x0没有风险约束，亦即风险率为0，因此使用以风险率为对角元素的矩阵来作为参数矩阵
A_eq = np.array([1.0])
A_eq = np.append(A_eq, 1 + df['p_i'] * 0.01)  # 此时只有一个等式约束，而A_eq要求传入二维矩阵
A_eq = np.array([A_eq, np.array([0.0, 0.0, 0.0, 0.0, 0.0])])  # 因此额外新增一行全零来凑够格式要求
b_eq = np.array([M, 0.0])  # 常数矩阵同样需要至少为两分量向量，令第二个分量为0，即可使第二个凑数约束恒成立
# 以上矩阵元素与a无关，不在循环时反复求取
for a in np.arange(0, 0.05, 0.001):  # 一系列a的值。从0开始，到0.05结束，步长为0.001
    b_ub = np.ones(5) * a * M  # 每个资产的风险都不能超过a*M，因此构造全1向量，再乘a*M
    b_ub[0] = 0.0  # x0代表银行里存的钱，是没有风险的，对其格外置一
    # 2.同样，如涉及不等式约束，也是需要让Ax<=b（左边为有变量的一边），需要对系数进行相同的处理。此处刚好符合形式，就不调整了
    # ————————————————————————————————调整————————————————————————————————
    result = linprog(c, A_ub, b_ub, A_eq, b_eq, [(0, None), (0, None), (0, None), (0, None), (0, None)])
    earnings = np.append(earnings, -1.0 * result.fun)
    earnings_pure = np.append(earnings_pure, -1.0 * result.fun - a * M)
    risk_thresholds = np.append(risk_thresholds, a)
    solutions.append(result.x)
plt.scatter(risk_thresholds, earnings)  # 画出收益随风险阈值变化的图像。此处观察原始图像后发现在较大时函数增率不明显，因此本来是想对它取对数的
# log()函数默认取自然对数为底，如果想取任意底数可以使用换底公式
# 读了值之后发现不是不明显，而是根本没在动了，哎算了不管它
plt.title("Risk thresholds && Earnings")
plt.xlabel("Risk thresholds")
plt.ylabel("Earnings")
plt.savefig('C:/Users/admin/Desktop/图一.png', dpi=500, bbox_inches='tight')
plt.show()
# 从图形可以看出来，从阈值略大于0.06开始，收益增长放缓，因此可以在从0到这个阈值的区间内进行选择
plt.scatter(risk_thresholds, earnings_pure)  # 然后是净利润随风险阈值变化的图像
plt.title("Risk thresholds && Pure earnings(the natural logarithm taken)")
plt.xlabel("Risk thresholds")
plt.ylabel("Earnings")
plt.savefig('C:/Users/admin/Desktop/图二.png', dpi=500, bbox_inches='tight')
plt.show()
index = np.argmax(earnings)
print(f"最高收益为{earnings[index]},对应风险阈值为{risk_thresholds[index]},（包括银行存款在内的）每种资产应分别投资{solutions[index]}")
print(f"转折点收益为{earnings[6]},对应风险阈值为{risk_thresholds[6]},（包括银行存款在内的）每种资产应分别投资{solutions[6]}")
# 如果想追求最大收益，定0.025的阈值即可，如果还想追求更加稳妥，就放低阈值。
# 以上的结论是十分主观的，如果想要达到理论、理性上的收益与风险的综合考量，可以采用TOPSIS法
# 先将风险转成最大型指标
max = np.max(risk_thresholds)
risk_thresholds = max - risk_thresholds


# 将两个指标进行标准化
def normalizing(data):
    return data / (data.dot(data) ** 0.5)


earnings = normalizing(earnings)
risk_thresholds = normalizing(risk_thresholds)

# 将两个向量放在一个矩阵里，方便进行后续处理
data = np.array([earnings, risk_thresholds]).T
max = np.max(data, axis=0)  # 取出每一行的最值，作为最优与最不优解，计算每个数据与它们的“距离”
min = np.min(data, axis=0)
result = np.array([])


def fun(data, max, min):
    d_pos = (data - max).dot(data - max) ** 0.5
    d_neg = (data - min).dot(data - min) ** 0.5
    return d_neg / (d_pos + d_neg)


data_len = len(data)
for i in range(data_len):
    result = np.append(result, fun(data[i], max, min))
print(f"所有解的得分为：{result},其中得分最高的结果对应风险阈值为{(np.argmax(result) + 1) * 0.001}")
# 阈值为0.07。看来选取拐点不完全是最优
